#This script executes regressions to estimate 
#associations between first-generation prenatal exposure to nonattainment 
#and first-generation time use.

rm(list=ls())
gc()
library(reldist)
library(ineq)
library(stargazer)
library(data.table)
library(dplyr)
library(dtplyr)
library(ggplot2)
library(stringdist)
library(stringr)
library(lfe)
library(haven)
library(broom)
library(tidylog)
library(foreign)
library(readr)

# Set the working directory
setwd("/projects/opp_env/caa1970/")

# Source external R scripts for custom functions
source("Code/rounding.r")

# Load the ATUS dataset

load("Data/ATUS/atus_merged_data.rda")

#Table 6

  Table_6 <- data.frame()
    
# Panel A: IV 

  Tab6_A1 <- felm(quality_time ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                    PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                    county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                    county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm|state_year+county +year_factor + month.y + survey_year_by_year + interview_day
                    |(all_fullterm~nonattainment_infant)|county, 
                  data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24, AGE < 46))
  
  summary(Tab6_A1)
  Tab6_A1$N
  
  temp_n<-  rounding_rules_counts(Tab6_A1$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24, AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time, na.rm=T ),4)
  temp_f <-  signif(Tab6_A1$stage1$iv1fstat[[1]]["F"],4)
  temp_results <- tidy(Tab6_A1) %>% filter(term == "`all_fullterm(fit)`" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl, fstat = temp_f,  table = "6", column ="A1" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  Tab6_A2 <- felm(quality_time_A ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                    PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                    county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                    county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm|state_year+county +year_factor + month.y + survey_year_by_year
                  |(all_fullterm~nonattainment_infant)|county, 
                  data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24, AGE < 46))
  
  summary(Tab6_A2)
  Tab6_A2$N
  temp_n<-  rounding_rules_counts(Tab6_A2$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24, AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time_A, na.rm=T ),4)
  temp_f <-  signif(Tab6_A2$stage1$iv1fstat[[1]]["F"],4)
  temp_results <- tidy(Tab6_A2) %>% filter(term == "`all_fullterm(fit)`" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl, fstat = temp_f,  table = "6", column ="A2" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  Tab6_A3 <- felm(quality_time_B ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                    PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                    county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                    county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm|state_year+county +year_factor + month.y + survey_year_by_year
                  |(all_fullterm~nonattainment_infant)|county, 
                  data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46))
  
  summary(Tab6_A3)
  Tab6_A3$N
  temp_n<-  rounding_rules_counts(Tab6_A3$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time_B, na.rm=T ),4)
  temp_f <-  signif(Tab6_A3$stage1$iv1fstat[[1]]["F"],4)
  temp_results <- tidy(Tab6_A3) %>% filter(term == "`all_fullterm(fit)`" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl, fstat = temp_f,  table = "6", column ="A3" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  
  
  Tab6_A4 <- felm(quality_time_C ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                    PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                    county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                    county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm|state_year+county +year_factor + month.y + survey_year_by_year
                  |(all_fullterm~nonattainment_infant)|county, 
                  data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE >25 , AGE < 46))
  
  summary(Tab6_A4)
  Tab6_A4$N
  temp_n<-  rounding_rules_counts(Tab6_A4$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time_C, na.rm=T ),4)
  temp_f <-  signif(Tab6_A4$stage1$iv1fstat[[1]]["F"],4)
  temp_results <- tidy(Tab6_A4) %>% filter(term == "`all_fullterm(fit)`" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl, fstat = temp_f,  table = "6", column ="A4" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  Tab6_A5 <- felm(quality_time_D ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                    PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                    county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                    county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm|state_year+county +year_factor + month.y + survey_year_by_year
                  |(all_fullterm~nonattainment_infant)|county, 
                  data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46))
  
  summary(Tab6_A5)
  Tab6_A5$N
  temp_n<-  rounding_rules_counts(Tab6_A5$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time_D, na.rm=T ),4)
  temp_f <-  signif(Tab6_A5$stage1$iv1fstat[[1]]["F"],4)
  temp_results <- tidy(Tab6_A5) %>% filter(term == "`all_fullterm(fit)`" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl, fstat = temp_f,  table = "6", column ="A5" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  Tab6_B1 <- felm(quality_time ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                    PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                    county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                    county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm+nonattainment_infant|state_year+county +year_factor + month.y + survey_year_by_year
                    |0|county, 
                  data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm)))
    

# Panel B: Reduced Form 

  summary(Tab6_B1)
  Tab6_B1$N
  temp_n<-  rounding_rules_counts(Tab6_B1$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time, na.rm=T ),4)
  temp_results <- tidy(Tab6_B1) %>% filter(term == "nonattainment_infant" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl,  table = "6", column ="B1" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  Tab6_B2 <-  felm(quality_time_A ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                     PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                     county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                     county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm+nonattainment_infant|state_year+county +year_factor + month.y + survey_year_by_year
                   |0|county, 
                   data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm)))
  
  summary(Tab6_B2)
  Tab6_B2$N
  temp_n<-  rounding_rules_counts(Tab6_B2$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time_A, na.rm=T ),4)
  temp_results <- tidy(Tab6_B2) %>% filter(term == "nonattainment_infant" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl,  table = "6", column ="B2" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  Tab6_B3 <-  felm(quality_time_B ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                     PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                     county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                     county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm+nonattainment_infant|state_year+county +year_factor + month.y + survey_year_by_year
                   |0|county, 
                   data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm)))
  
  summary(Tab6_B3)
  Tab6_B3$N
  temp_n<-  rounding_rules_counts(Tab6_B3$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time_B, na.rm=T ),4)
  temp_results <- tidy(Tab6_B3) %>% filter(term == "nonattainment_infant" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl,  table = "6", column ="B3" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  Tab6_B4 <-  felm(quality_time_C ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                     PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                     county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                     county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm+nonattainment_infant|state_year+county +year_factor + month.y + survey_year_by_year
                   |0|county, 
                   data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm)))
  
  summary(Tab6_B4)
  Tab6_B4$N
  temp_n<-  rounding_rules_counts(Tab6_B4$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time_C, na.rm=T ),4)
  temp_results <- tidy(Tab6_B4) %>% filter(term == "nonattainment_infant" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl,  table = "6", column ="B4" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  
  Tab6_B5 <-  felm(quality_time_D ~factor(race.x)+factor(sex.x)+PI_PC*year_trend+
                     PI_PC*I(year_trend^2) +epop*year_trend + total_transfers*year_trend + 
                     county_pop*year_trend + epop*I(year_trend^2) + total_transfers*I(year_trend^2) + 
                     county_pop*I(year_trend^2) +tmax_fullterm+ tmin_fullterm+prec_fullterm+anyprec_fullterm+nonattainment_infant|state_year+county +year_factor + month.y + survey_year_by_year
                   |0|county, 
                   data=cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm)))
  
  summary(Tab6_B5)
  Tab6_B5$N
  temp_n<-  rounding_rules_counts(Tab6_B5$N)
  temp_ctrl <- signif(mean((cps_ipums_atus_merge_pob_caa %>% filter(year.y%in%1969:1980, AGE>24 , AGE < 46, !is.na(all_fullterm),nonattainment_infant == 0))$quality_time_D, na.rm=T ),4)
  temp_results <- tidy(Tab6_B5) %>% filter(term == "nonattainment_infant" ) %>% #tidy regression, keep coef of interest
    mutate_each(funs(signif(.,4)), -term) %>% #round appropriately
    mutate(N=temp_n, ctrl_mean=temp_ctrl,  table = "6", column ="B5" ) #add variables for N and control mean, plus index the table number and column
  Table_6 <- bind_rows(Table_6, temp_results) #Appended to big dataset
  write_csv(Table_6, "/projects/opp_env/caa1970/JPE_Micro/Output/Table 6/Table_6.csv")
  

